from abc import ABC, abstractmethod
import math
from collections import deque
import shapely
from ..data_structures.base_structures import *

class BaseQuantifier(ABC):
    """量化器基类"""
    @abstractmethod
    def quantify_hexcells(path:str, map:Map):
        """量化HexCell"""
        pass

    @abstractmethod
    def quantify_quadcells(path:str, map:Map):
        """量化QuadCell"""
        pass

    def geohash_polyfill_geojson(geojson, precision):
        """
        使用邻接递推方式填充 GeoJSON Polygon 区域，生成紧密、无重复的 geohash 编码集合。
        :param geojson: GeoJSON 格式（type 为 Polygon 或 Feature）
        :param precision: geohash 编码精度
        :return: Set[str] 被覆盖的 geohash 编码集合
        """
        def geohash_to_polygon(gh):
            """获取 geohash 对应的 Polygon 区域"""
            b = geohash.bbox(gh)
            return shapely.geometry.Polygon([
                (b['w'], b['s']), (b['e'], b['s']),
                (b['e'], b['n']), (b['w'], b['n']),
                (b['w'], b['s'])  # 封闭
            ])

        # 提取 Polygon 对象
        if geojson.get("type") == "Feature":
            coordinates = geojson["geometry"]["coordinates"]
        elif geojson.get("type") == "Polygon":
            coordinates = geojson["coordinates"]
        else:
            raise ValueError("geojson 输入必须是 Feature 或 Polygon")

        poly = shapely.geometry.Polygon(coordinates[0])

        # 获取 bbox 中心的 geohash 编码
        minx, miny, maxx, maxy = poly.bounds
        center_lat = (miny + maxy) / 2
        center_lon = (minx + maxx) / 2
        start_gh = geohash.encode(center_lat, center_lon, precision)

        # 递归填充 geohash 编码
        def find_valid_seed(start_gh, geometry, max_depth=50):
            """
            若中心 geohash 无交集，向邻居扩展，寻找第一个与 polygon 相交的 geohash。
            """
            visited = set()
            queue = deque([start_gh])
            depth = 0
            while queue and depth <= max_depth:
                next_queue = deque()
                while queue:
                    gh = queue.popleft()
                    if gh in visited:
                        continue
                    visited.add(gh)
                    gh_poly = geohash_to_polygon(gh)
                    if geometry.intersects(gh_poly):
                        return gh  # 找到第一个相交的 gh
                    next_queue.extend([n for n in geohash.neighbors(gh) if n not in visited])
                queue = next_queue
                depth += 1
            return None  # 未找到合适起点

        # 3. 检查中心是否有效，否则递归寻找邻居
        if not poly.intersects(geohash_to_polygon(start_gh)):
            start_gh = find_valid_seed(start_gh, poly)
            if start_gh is None:
                return set()

        # 4. BFS 方式向外填充 geohash
        visited = set()
        result = set()
        queue = deque([start_gh])

        while queue:
            gh = queue.popleft()
            if gh in visited:
                continue
            visited.add(gh)

            gh_poly = geohash_to_polygon(gh)
            if poly.intersects(gh_poly):
                result.add(gh)
                for ngh in geohash.neighbors(gh):
                    if ngh not in visited:
                        queue.append(ngh)

        return result

    def estimate_cell_size(precision):
        """
        根据 geohash 精度计算其对应的大致边长（单位：米）
        :param precision: geohash 精度（整数）
        :return: 边长（单位：米）
        """
        precision = int(precision)  # 强制转换为整数
        if not (1 <= precision <= 12):
            raise ValueError("Geohash 精度应为 1~12 之间的整数")

        enum_name = f"PRECISION_{precision}"
        width, height = GeohashCellSize[enum_name].value
        # 计算平均边长
        cell_size = (width + height) / 2
        # 返回平均边长（也可以返回宽高元组）
        return cell_size

    def lonlat_to_pixel(transform, lon, lat):
        """将经纬度转换为栅格的行列索引"""
        col, row = ~transform * (lon, lat)
        return int(row), int(col)

    def get_value(data, transform, lon, lat, nodata_value):
        """根据经纬度获取高程"""
        row, col = BaseQuantifier.lonlat_to_pixel(transform, lon, lat)
        if 0 <= row < data.shape[0] and 0 <= col < data.shape[1]:
            value = float(data[row, col])
            return value if value != nodata_value else None  # 排除缺省值
        return None  # 超出边界返回 None
    
    def is_within_bounds(lat, lon, bounds):
        """判断经纬度是否在栅格范围内"""
        min_lon, min_lat = bounds.left, bounds.bottom
        max_lon, max_lat = bounds.right, bounds.top
        return min_lon <= lon <= max_lon and min_lat <= lat <= max_lat
    
    def get_neighbor_elevations(map:Map, cell):
        """
        获取cell及其邻居的高程值,组成3x3网格
        兼容四角格网（8邻）与六边形格网（6邻）
        :param map: 地图对象
        :param cell: 当前cell
        :return: 3x3高程数组
        """
        # 获取当前cell的高程
        center_elevation = cell.get_attribute(DemAttribute.Elevation).get_value()
        
        # 初始化3x3网格，中心为当前cell
        elevations = [None] * 9
        elevations[4] = center_elevation  # 中心位置
        
        # 获取邻居cell的高程
        neighbors = cell.neighbors
        if len(neighbors) == 6:  # H3六边形有6个邻居
            # 将6个邻居映射到3x3网格的8个位置
            # 这里使用简化映射，将邻居按照相对位置分布
            for i, neighbor_idx in enumerate(neighbors):
                neighbor_cell = map.cells.get(neighbor_idx)
                if neighbor_cell is None:
                    continue
                neighbor_cell_elevation = neighbor_cell.get_attribute(DemAttribute.Elevation).get_value()
                if neighbor_cell and neighbor_cell_elevation is not None:
                    # 将6个邻居映射到3x3网格的相应位置
                    grid_pos = [0, 1, 2, 5, 8, 7, 6, 3][i] if i < 8 else 4
                    elevations[grid_pos] = neighbor_cell_elevation
        if len(cell.neighbors) == 8:  # 四角格网有8个邻居
            neighbor_mapping = [0, 1, 2, 5, 8, 7, 6, 3]  # 顺时针
            # 遍历邻居编码并填充高程值
            for i, neighbor_idx in enumerate(cell.neighbors[:8]):
                neighbor_cell = map.cells.get(neighbor_idx)
                if neighbor_cell == None:
                    continue
                neighbor_cell_elevation =neighbor_cell.get_attribute(DemAttribute.Elevation).get_value()
                if neighbor_cell and neighbor_cell_elevation is not None:
                    elevations[neighbor_mapping[i]] = neighbor_cell_elevation
        
        return elevations
    
    def calculate_curvature(elevations, cell_size):
        """
        计算曲率(使用Horn算法的二阶导数)
        :param elevations: 3x3高程数组
        :param cell_size: 单元格大小(米)
        :return: 曲率值
        """
        # 3x3窗口布局:
        # a b c
        # d e f  
        # g h i
        a, b, c, d, e, f, g, h, i = elevations
        
        # 如果中心点为None，返回None
        if e is None:
            return None
            
        # 处理缺失值，用中心值替代
        def fill_none(val, default):
            return default if val is None else val
            
        a = fill_none(a, e)
        b = fill_none(b, e)
        c = fill_none(c, e)
        d = fill_none(d, e)
        f = fill_none(f, e)
        g = fill_none(g, e)
        h = fill_none(h, e)
        i = fill_none(i, e)
        
        # 计算一阶偏导数
        dz_dx = ((c + 2*f + i) - (a + 2*d + g)) / (8 * cell_size)
        dz_dy = ((g + 2*h + i) - (a + 2*b + c)) / (8 * cell_size)
        
        # 计算二阶偏导数
        d2z_dx2 = ((a + 2*b + c) + (g + 2*h + i) - 2*(d + 2*e + f)) / (3 * cell_size * cell_size)
        d2z_dy2 = ((a + 2*d + g) + (c + 2*f + i) - 2*(b + 2*e + h)) / (3 * cell_size * cell_size)
        d2z_dxdy = ((a + i) - (c + g)) / (4 * cell_size * cell_size)
        
        # 计算平均曲率 (Mean Curvature)
        p = dz_dx
        q = dz_dy
        r = d2z_dx2
        s = d2z_dxdy
        t = d2z_dy2
        
        # 避免除零错误
        denominator = (1 + p*p + q*q)
        if denominator == 0:
            return 0
            
        # 平均曲率公式
        mean_curvature = -((1 + q*q) * r - 2*p*q*s + (1 + p*p) * t) / (2 * (denominator ** 1.5))
        
        return round(mean_curvature * 1000, 6)  # 放大1000倍便于观察，保留6位小数

    def calculate_gradient(elevations, cell_size):
        """
        计算3x3窗口的梯度
        :param elevations: 3x3高程数组
        :param cell_size: 单元格大小(米)
        :return: (dz_dx, dz_dy) 梯度值
        """
        # Horn算法计算梯度
        # 3x3窗口布局:
        # a b c
        # d e f  
        # g h i
        a, b, c, d, e, f, g, h, i = elevations
        
        # 如果中心点或关键点为None，返回None
        if e is None:
            return None, None
            
        # 处理缺失值，用中心值或邻近值替代
        def fill_none(val, default):
            return default if val is None else val
            
        a = fill_none(a, e)
        b = fill_none(b, e)
        c = fill_none(c, e)
        d = fill_none(d, e)
        f = fill_none(f, e)
        g = fill_none(g, e)
        h = fill_none(h, e)
        i = fill_none(i, e)
        
        # Horn算法
        dz_dx = ((c + 2*f + i) - (a + 2*d + g)) / (8 * cell_size)
        dz_dy = ((g + 2*h + i) - (a + 2*b + c)) / (8 * cell_size)
        
        return dz_dx, dz_dy

    def calculate_exposure_from_gradient(dz_dx, dz_dy):
        """
        根据梯度计算坡向
        :param dz_dx: x方向梯度
        :param dz_dy: y方向梯度
        :return: 坡向角度(度),北为0度,顺时针增加
        """
        if dz_dx == 0 and dz_dy == 0:
            return None  # 平地无坡向
        
        # 计算坡向角度（弧度）
        aspect_rad = math.atan2(dz_dy, -dz_dx)
        
        # 转换为度数，并调整为北为0度的坐标系
        aspect_deg = math.degrees(aspect_rad)
        aspect_deg = 90 - aspect_deg
        
        # 确保角度在0-360度范围内
        if aspect_deg < 0:
            aspect_deg += 360
        elif aspect_deg >= 360:
            aspect_deg -= 360
            
        return round(aspect_deg, 2)
    
    def calculate_center_slope(hexcell_size, vertex_elevation, center_elevation, grid_type='hex'):
        """使用Horn算法计算H3六边形单元的坡度"""
        if center_elevation is None :
            return None

        # 判断顶点数量是否匹配
        if grid_type == 'hex' and len(vertex_elevation) != 6:
            return None
        elif grid_type == 'quad' and len(vertex_elevation) != 4:
            return None

        # 构建3x3高程矩阵，中心为当前单元格中心点
        elev_matrix = [[None, None, None],
                      [None, center_elevation, None],
                      [None, None, None]]
                      
        # H3六边形顶点在3x3矩阵中的位置映射（顺时针排列，从顶部开始）
        # 0: 上, 1: 右上, 2: 右下, 3: 下, 4: 左下, 5: 左上
        if grid_type == 'hex':
            vertex_positions = [
                (0, 1),  # 上
                (0, 2),  # 右上
                (2, 2),  # 右下
                (2, 1),  # 下
                (2, 0),  # 左下
                (0, 0)   # 左上
            ]
        else:  # 四角格网
            vertex_positions = [
                (0, 0),  # 左上
                (0, 2),  # 右上
                (2, 2),  # 右下
                (2, 0),  # 左下
            ]
        
        # 填充顶点高程值
        valid_points = 0
        for i, elev in enumerate(vertex_elevation):
            if elev is not None:
                row, col = vertex_positions[i]
                elev_matrix[row][col] = elev
                valid_points += 1
                
        # 需要至少4个有效点（包括中心点）才能计算坡度
        if valid_points < 3:  # 3个顶点 + 1个中心点
            return None
            
        # 计算x方向和y方向的梯度
        # Horn算法的权重系数
        dz_dx = (
            (elev_matrix[0][2] if elev_matrix[0][2] is not None else center_elevation) +
            2 * (elev_matrix[1][2] if elev_matrix[1][2] is not None else center_elevation) +
            (elev_matrix[2][2] if elev_matrix[2][2] is not None else center_elevation) -
            (elev_matrix[0][0] if elev_matrix[0][0] is not None else center_elevation) -
            2 * (elev_matrix[1][0] if elev_matrix[1][0] is not None else center_elevation) -
            (elev_matrix[2][0] if elev_matrix[2][0] is not None else center_elevation)
        ) / (8.0 * hexcell_size)
        
        dz_dy = (
            (elev_matrix[2][0] if elev_matrix[2][0] is not None else center_elevation) +
            2 * (elev_matrix[2][1] if elev_matrix[2][1] is not None else center_elevation) +
            (elev_matrix[2][2] if elev_matrix[2][2] is not None else center_elevation) -
            (elev_matrix[0][0] if elev_matrix[0][0] is not None else center_elevation) -
            2 * (elev_matrix[0][1] if elev_matrix[0][1] is not None else center_elevation) -
            (elev_matrix[0][2] if elev_matrix[0][2] is not None else center_elevation)
        ) / (8.0 * hexcell_size)
        
        # 计算坡度（度）
        slope_rad = math.atan(math.sqrt(dz_dx * dz_dx + dz_dy * dz_dy))
        slope_deg = math.degrees(slope_rad)
        
        return slope_deg


